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Abstract 

We investigate the out of equihbrium dynamics of the two-dimensional 
XY model when cooled across the Berezinskii-Kosterlitz-Thouless (BKT) 
phase transition using different protocols. We focus on the evolution of the 
growing correlation length and the density of topological defects (vortices) . 
By using Monte Carlo simulations we first determine the time and tem- 
perature dependence of the growing correlation length after an infinitely 
rapid quench from above the transition temperature to the quasi-long 
range order region. The functional form is consistent with a logarithmic 
correction to the diffusive law and it serves to validate dynamic scaling in 
this problem. This analysis clarifies the different dynamic roles played by 
bound and free vortices. We then revisit the Kibble-Zurek mechanism in 
thermal phase transitions in which the disordered state is plagued with 
topological defects. We provide a theory of quenching rate dependence in 
systems with the BKT-type transition that goes beyond the equilibrium 
scaling arguments. Finally, we discuss the implications of our results to 
a host of physical systems with vortex excitations including planar ferro- 
magnets and liquid crystals as well as the Ginzburg-Landau approach to 
bidimensional freely decaying turbulence. 

1 Introduction 

The out of equilibrium dynamics of systems with continuous symmetries an- 
nealed or quenched from the symmetric to the symmetry broken phase are 
a subject of study in different branches of physics. Liquid-crystals quenched 
across a phase transition into their ordered phase present a large variety of 
topological defects that diffuse, interact and eventually annihilate [1 . Vector 
ferromagents cooled across their Curie temperature are other condensed-matter 
examples with a variety of dynamic topological defects [2 . Cosmology provides 
a set of interesting models in which similar dynamics occur [3 . Still, the decay 



of vortex density in freely decaying turbulence has also been of interest [T. In 
all these areas, there is considerable interest in the dynamics of defects once the 
phase transition has been crossed. 

Of special interest is the planar time-dependent Ginzburg Landau model 
with SO (2) symmetry or its lattice counterpart the 2d XY model. A static 
phase transition occurs at a finite critical temperature, T^t, between a high-T 
disordered paramagnet and a low-T ferromagnetic phase with quasi-long range 
order. When the model is quenched from the high to the low temperature 
phase the continuous symmetry of the disordered phase is broken to the one 
of the ordered phase. In these systems the defects are singular vortices that 
carry topological charge and have logarithmic interactions. In the high-T phase 
free vortices proliferate while in the low-T phase vortices bind in pairs. Physical 
realizations are bidimensional planar ferromagnets [2] , superconducting films [5] , 
Josephson-j unction arrays [6 , especially tailored nematic liquid crystals [7], and 
toy models for bidimensional turbulence [4]. The exponential singularity may 
also describe the critical properties of superfiuid helium films [8 . 

Phase ordering kinetics following an infinitely rapid quench through a ther- 
mal critical point have been extensively studied in the statistical physics liter- 
ature [9 . In condensed matter applications the time spent close to the critical 
point is typically much shorter than the time spent far from it and, therefore, 
the infinitely fast quench approximation is justified. A few studies where cooling 
rate dependencies have been taken into account are [lOl [11] [12] [13] . The study 
of cooling-rate dependencies in defect dynamics is, instead, central in cosmol- 
ogy. The goal in this case is to estimate the density of defects, p, left over in 
the universe after it went through a number of potential (second-order) phase 
transitions. The importance lies on the fact that, initially, these were thought 
to act as seeds for matter clustering. This scenario seems to be excluded by 
observation data now but, still, the interest in predicting the density of topo- 
logical defects remains due to other possible effects of these objects [14^. In the 
late 80s Zurek derived a quantitative prediction for p that is based on critical 
scaling in the disordered phase close to a second order phase transition. The 
proposal is based upon the hypothesis whereby the defect dynamics should be 
negligible below the transition. In consequence, p is assumed to remain fixed to 
the value it takes at the control parameter at which the system falls out of equi- 
librium in the disordered phase. In his articles Zurek proposed to check these 
predictions in condensed matter systems with the same symmetry properties as 
the cosmological models [15 . A vast experimental [161 [l^ numerical 
activity followed with variable results summarized in fi9\. 

The fact that the annihilation of defects can be neglected in the out of 
equilibrium regime was questioned in [20 . In this paper, a scaling argument 
that includes the mechanism of defect annihilation in the low-temperature phase 
for systems with second-order phase transitions and dissipative dynamics was 
proposed. A scaling of the density of defects with cooling rate and time spent 
in the out of equilibrium evolution was derived. The predictions were checked 
numerically with Monte Carlo dynamics of the 2d Ising model, a paradigmatic 
system with discrete symmetry breaking. Similar ideas about the importance 
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of the dynamics below the transition were stressed in [21]. 

The treatment of ah cases mentioned in the previous paragraph is classical. 
Recently, pushed by the advent of powerful experimental techniques in cold 
atom systems, the extension of the Kibble-Zurek mechanism to quantum isolated 
systems was developed [22j [23l [24] . These claims have been critically revisited 
in [25l [26] , as reviewed in [27] . 

In this paper we examine the density of vortices left over after going through 
the Berezinskii-Kosterlitz-Thouless classical phase transition ^28^ ^29^ with a fi- 
nite speed. We revisit the Kibble-Zurek mechanism and we provide a theory 
of quenching rate dependence in systems with continuous symmetry and the 
BKT-type of 'infinite-order' phase transition that goes beyond the equilibrium 
scaling arguments. We test it with a Monte Carlo numerical study of the 2d 
XY model on a square lattice. We work in the canonical setting, in the sense 
that the system is coupled to an equilibrium environment that allows for energy 
dissipation. As a control parameter driving the phase transition we use the tem- 
perature of the heat-bath. For concreteness we use linear cooling procedures, 
characterized by a single parameter, the cooling rate. A previous numerical 
study of this problem with a different method that ignores the effects of free- 
vortices appeared in [30^ . We shall discuss the results in this paper and confront 
them to ours in the body of this work. 

The organization of the paper is the following. In Sect. |2] we recall the defi- 
nitions of the 2d XY model and the time-dependent Ginzburg-Landau equation 
and we present the Monte Carlo method used. In Sect. |3] we shortly review the 
BKT scenario. Section |4] is devoted to the analysis of infinitely rapid quenches. 
In Sect. |5] we present our results on annealing procedures. Finally, in Sect. [6] we 
give our conclusions. 

2 Model and method 

In this section we recall the definition of the 2d XY model and we describe 
the Monte Carlo method that we use in our simulations of the lattice model. 
We also explain the corresponding time-dependent Ginzburg-Landau stochastic 
equation. 

2.1 The 2d XY model 

The two-dimensional XY model is defined by the Hamiltonian 

H=-JY,S^■S,, (1) 

with ferromagnetic exchange coupling, J > 0. The spins are classical variables 
constrained to live on a plane, sf = 1, and the sum runs over nearest neighbours 
on a square lattice of linear size L. It is convenient to use a parametrization in 
which the spin vector is represented by the angle it forms with a chosen axis. 
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= ; with it the Hamiltonian becomes 

H = - J^cos((9i -6>^) . (2) 

Above the critical temperature T^t, unbound positively and negatively charged 
vortices are present and the system is disordered. Below Tkt vortex- antivortex 
pairs bind and the system is critical in a sense that we shall make more precise 
below. The critical temperature is found to be ksTxT — 0.89 J [29 . We measure 
T in units of J /ks in the following. 

2.2 Dynamics 

Stochastic dynamics, mimicking the coupling of the classical spins to a thermal 
bath, can be attributed in different ways. In this Section we define the Monte 
Carlo procedure used in our simulations. We focus on non-conserved order- 
parameter dynamics (model A) in which the magnetization is not conserved. 
We also briefiy recall the time-dependent Ginzburg-Landau approach. 

2.2.1 Monte Carlo dynamics 

We performed a Monte Carlo study of the 2d XY model on a square lattice of 
linear size L with periodic boundary conditions. The Monte Carlo algorithm 
consists in updating the angular variable Oi associated to a randomly chosen 
site z to a new value 9[ randomly chosen in the interval [— 7r,7r], with proba- 
bility p = min {l, e~^^^}, where /S.E is the energy variation between the two 
configurations. Equilibrium data are for relatively small system sizes, L = 16 
to L = 50, while out of equilibrium ones are for larger samples, L = 100 to 
L = 400. We briefiy discuss finite size effects when presenting the simulation 
data. The results are an average of 5 to 10000 independent runs each (depend- 
ing on the quantity calculated) and one unit of time is an attempted move of 
every spin. 

In infinitely rapid quenches from infinite temperature we took a fully random 
initial condition in which a random number Oi in the interval [— tt, tt] is asso- 
ciated to site i. We explain the implementation of more sophisticated cooling 
procedures in Sect.jS] 

The focus of our study is the evolution of the number of vortices. We de- 
termine it numerically with two types of measurements. On the one hand dy- 
namic scaling implies that the density of vortices, p^, should depend on the 
tyical growing lengh, ^, as py We shall use several determinations 

of ^ that include the exponential decay of the two-point correlation function, 
C[(,c{t^T)^t] = const., and the second moment of C{r,t), that is defined as 
follows [31 . From the total magnetization M{t) = Xl^Si(t) and the magnetic 
susceptibility 

X{t) = ^M'{t), (3) 
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the second moment correlation length on a lattice of size is defined by 
where 

^(*) = ^^'^^^ ' ^'^^^^ ^^''^^^ " '"^^ ■ 

In addition, we determine the number of vortices Ny directly from the con- 
figuration of the system by counting the plaquettes with non-zero vorticity. The 
vorticity is determined as the integer winding number n, so that the phase dif- 
ference Oij = Oi — Oj around the plaquette is equal to ^ Oij = 27rn. The number 
of vortices Ny is then obtained as the number of plaquettes around which the 
phase 6i rotates through ±27r, taking care that the phase difference % is re- 
stricted to the interval [— 7r,7r], as are the phases Oi. The density of vortices is 
then py = Ny/L'^. 

In the figures we will not report the statistical errors because they are very 
small for the quantities we consider (smaller than dot dimensions in the vari- 
ous plots). For example, for the largest system size (L = 400) for which the 
fluctuations are larger, the variance of the number of vortices is at most 3 to 
4%. 



2.2.2 Time-dependent Ginzburg-Landau equation 

A field theory version of this problem, better suited for analytic calculations 
is given by the time-dependent Ginzburg-Landau equation that determines the 
stochastic evolution of the coarse-grained two-component field (j): 

dkr,t) ^ ^ SH 
dt 5${r,t) 
with r the kinetic coefficient. The 'free-energy' is 



H 



j d'r[e^(Smf + \{4>\^-lf] (7) 



with ps the stiffness coefficient or elastic constant and u> {) the strength of the 
non-linearity that drives the vector field to take unit modulus. Noise can be 
added to Eq. ([6| in the form of an additional stochastic term in its right-hand- 
side with Gaussian statistics and, typically, delta-correlated, {Ca{^^t)Cb{^ ^t')) = 
2kBTT5ah^if—f')5{t — t')^ with a, 6 = 1, 2. The analytic predictions on the time- 
dependent growing length to be discussed in Sect. |4.l| have been obtained using 
this approach [32j [33] . 



3 The BKT equilibrium phase transition 

We summarize the static properties of the 2d XY model as figured out by 
Berezinskii [28 and Kosterlitz and Thouless [29 (the BKT picture). 
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In the thermodynamic hmit 2d XY models display no conventional long- 
range order (the magnetization vanishes at all temperatures due to spin-wave 
excitations) but topological order below a critical temperature Tkt- The physics 
is dominated by two types of excitations: harmonic spin waves and vortices. 
The former are responsible for destroying conventional long-range order while 
the latter are responsible for the phase transition. The Berezinskii-Kosterlitz- 
Thouless (BKT) phase transition at Tkt is characterized by a change in the 
behaviour of the equilibrium spatial correlation function C{r) = _^ 
where the angular brackets denote an average over the equilibrium measure. 
The high-T phase is disordered, the correlation decays exponentially, 

C{r) - e-^/^^^ , (8) 

and there is a finite density of free vortices. Close and above Tkt the correlation 
length diverges exponentially 

Ce^^aj e^sK^-^^^)/^'^^!"" (9) 

with u = 1/2 and a non- universal constant that typically takes a value of 
order one (on a square lattice ~ 1.5) [28l[29]. Thermodynamic quantities are 
smooth accross the transition. The low-T phase has quasi-long-range order, 

C{r) ^ r-^^^-^m = r-^m , (10) 

it is critical in the sense that the equilibrium correlation length diverges, 
oo, and it is characterized by the T-dependent exponent r]{T) [29] that decreases 
upon decreasing temperature from r]{TKT) = 1/4 to r]{T ^ 0) = T/(27rJ). 
Bound vort ex- ant i vortex pairs populate the low-temperature ordered phase co- 
exiting with spin- waves. The transition is interpreted by using an analogy with 
the Coulomb gas in which the system has a low-temperature dielectric phase 
(with charges, in the 2d XY case vortices, bound into dipoles) and a high- 
temperature plasma or conducting phase with free charges (in the 2d XY model 
vortices) [28", "29^. In short, at Tkt pairs dissociate. 

The BKT singularity yields the best fit to data generated with large-scale 
Monte Carlo simulations and short-time dynamics. The data tend to confirm 
the analytical values of the exponents u and r]{T) and the critical temperature 
was estimated to Tkt — 0.89 [34J [SS] |36l |371 138]. In finite size systems the 
effective transition temperature behaves as Tkt + 0(ln~^ L) |39]. Monte Carlo 
values of the correlation length £,eq{T) obtained with lattices with linear size 
L = 512 are Ceg(1.25) 3.8, ^eg(1.04) ^ 18.7, and ^eg(0.98) 70 We shah 



discuss our own data, obtained with much smaller system sizes, in Fig. 11 



4 Relaxation after a quench 

The nonequilibrium behavior of the 2d XY model following an instantaneous 
quench below Tkt has been studied theoretically as well as experimentally [71 



6 



|4Q] . In this Section we recall the time-dependence of the growing correlation- 
length. We go beyond this by-now well-established aspect of the dynamics with 
a careful analysis of the temperature-dependence that will be of use in the rest 
of the article. We present a detailed analysis of the vortex- ant i- vortex time- 
dependent structure. 



4.1 Analytic prediction for the growing length 

In this Section we set the origin of time to be the instant at which the initial 
condition is let evolve in the new parameter conditions. At all T < Tkt a system 
prepared in an out of equilibrium initial state approaches equilibrium through 
a coarsening process in which local - critical - equilibrium is established over 
a length scale ^(t,T). The space-time correlation C(r, t) = {si{t)Sj{t))\^^ _^^^^^ 
is governed by the scaling form 

= (^)^ (11) 

ri(T) is the static critical exponent. The scaling function depends on the initial 
conditions. We focus on high-temperature, disordered ones. /(O) = 1 so that 



the equilibrium result (10) is recovered for all a <C r <C ^(t,T), a being a 
microscopic cut-off length, and, in particular, in the infinite long time limit 
in which ^(t,T) oo. In the opposite limit r ^ ^{t^T) the system remains 
disordered, as in the initial condition, and this is ensured by an f{x) that rapidly 
falls-off for X > 1. 

The relaxation of an out of equilibrium state is due to two processes: the 
dynamics of Goldstone modes and the annihilation and pairing of vortices. One 
can tune the importance of the vortex contribution by choosing the initial con- 
dition. For example, a completely ordered configuration is free of vortices while 
a high temperature configuration has a finite density of free vortices. After 
taking a fully ordered initial condition to < T < Tkt^ no free vortices are 
generated by thermal fiuctuations and the system orders in the new conditions 
by growing a length ^(t,T) cx [Ao(T)t]^/^ [411 [42]. Instead, the relaxation of an 
initial condition with free vortices in the critical low-temperature phase is far 
more interesting. Whether dynamic scaling holds and, in the affirmative, which 
is the time-dependent growth-law was a subject of debate for some time. The 
following scenario is well-established now. 

The numerical integration of the zero-temperature Ginzburg-Landau equa- 
tion using fully random initial conditions (with, on average, a density of vortices 
of 1/3) first suggested a power-law for ^(t,T) with a small but measurable de- 
viation from the diffusive law t^/^ [43l[44]. This regime was considered to be a 
transient fully determined by the vortex annihilation dynamics. Yurke et al. ^] 
and Bray et al. [33 built upon these results and gave an argument for 



l/z 



(12) 
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in quenches from equilibrium at T > Tkt to T = 0, see below. For the sake of 
completeness, we wrote this expression inserting the dynamic exponents z {= 2 
in this case). The extension to finite working T was worked out in [45]. This law 
was confirmed by extensive zero-temperature numerical simulations in [46l |47] . 
The parameter A(T) and the microscopic cut-off to(^) are non- universal in the 
sense that they depend on the procedure used to measure ^. 



Let us briefly reproduce here the argument that leads to ( 12 ) and extend it 
later to infer the relaxation time. Within the Ginzburg-Landau model ([6|, a 
field configuration, describing a single free vortex (j) = f/r has an energy Ey = 
7rpsln(L/a), where L and a are the system size and microscopic cut-off (e.^., 
the lattice spacing). In the many- vortex situation, where ^{t) represents the 
typical spacing between vortices and anti-vortices one can derive the following 
expression for the typical speed of a vortex [32] [45] 

dm ^ p,r 

from the zero temperature time-dependent Ginzburg-Landau equation (|6|. The 



solution is Eq. \12\ with to — o? /ps^- We assume, following [45 , that this 
equation holds at finite T, and even above Tkt where, though, the growing 
length saturates to the equilibrium correlation length at sufficiently long times. 
Matching the growing length with the equilibrium law by d^/dt\^^ £,eq/'Teq 
above Tkt one finds 

re,~C:,ln(W«) • (14) 

An exponential divergence of Tgg with the distance to criticality is recovered by 
substituting with the expression in (|9|: 

re, ^ a. e^c^[(^-^K-)/T..]-^ f . (15) 

V -l-KT J 

The same behavior is estimated in [48^ with a different argument. 



4.2 Numerical measurements of the growing length 

Although different determinations of the growing length should converge, asymp- 
totically to the same and unique law, some are more convenient than others at 
finite times. Comparison of different prescriptions (from scaling of the space- 
time correlation, the density of topological defects, and the structure factor) 
were discussed in [46l [47]. We revisit this issue with the purpose of clarify- 
ing the different role played by bound and free vortices and establishing the 
temperature dependence of the correlation length. 

In Fig. [l] we extract the growing length ^(t,T) from the decay of the space- 
time correlation. In the left panel we show the bare data, and we confirm 
that C satisfies dynamic scaling by using the value of ^c{t^T) inferred from 
the condition C[£,c{t^T)^t] = 0.3. In the right panel we compare the time- 
dependence of ^c{t^T) to the power law t^/^ and we see the expected slower 
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Figure 1: (Color online.) Space-time correlation after an infinitely rapid quench 
from a fully random initial condition into the low temperature phase, (a) Bare 
data for C(r, t) at different times t after a quench to T = 0.5 that are given in the 
key. The inset shows the scaling plot r^^^^C against r/^^dt^T) using ^c{tjT) 
computed from C[£,c{t^T)^t] = 0.3. (b) ^c{t,T) against t for instantaneous 
quenches to different temperatures. 



growth at sufficiently long times. The small deviations due to the logarithmic 
correction will be studied in more detail in Fig. [2] 

Next we analyze the growing length as extracted from the density of vor- 
tices. In the fully random, infinite temperature limit, :^ 0.3. In equilibrium 
at 2Tkt we found a value slightly smaller than 0.2 in agreement with data shown 
in [35]. After the infinitely rapid quench the very high density of free vortices 
present in the initial high temperature condition decreases. Two processes con- 
tribute to the re-organization of the vortex configuration. On the one hand 
over-abundant defects annihilate. On the other hand vortices and anti-vortices 
bind to approach the non-vanishing equihbrium density of pairs, peq — e~^^/^ 
with II 3.77 [35], that is an increasing function of T. 

According to dynamic scaling, during the out of equilibrium relaxation the 
density of defects in the ordered phase should decrease as p^(t,T) 1/^^(LT) 
that using eq. (12) is equivalent to py{t^T)\n py{t^T) [32^. In Fig. 0we 

show our data for all vortices in two ways: in the left panel we plot pv{t) as a 
function of t for several values of T; in the right panel Ny (t) In Ny (t) as a function 
of t where Ny is the number of vortices selecting the case T = 0.4. After one 
MC step the density is very close to the initial value py{T oo) 0.3 but it 
subsequently monotonically decreases in time in all cases. As already stressed by 
Yurke et al [32] , who solved numerically the time-dependent Ginzburg-Landau 
equation with noise, the second presentation gives a much better description 
of data for T ~ 0.4, confirming the growth-law (12). However, at T close to 
Tkt the number of vortices gets close to the equilibrium value, for example 
Peq — 10~^ at T = 0.8, and the cross-over to equilibrium dynamics is reached 
within the simulation times. This is accompanied by the generation of many 
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short-lived vortex anti- vortex pairs, similar to what was observed in [32], and the 
data deviates from the expected scaling. The failure of scaling of the space-time 
correlation when is used was also mentioned in [46l |47S . 
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Figure 2: (Color online.) Time-dependence of the total density of vortices, 
pyit^T) = Nyit^T) / L'^ ^ after an infinitely rapid quench from Tq ^ oo to < 
T < Tkt given in the key. The data are shown in the form pv{t) vs. t in (a) and 
Ny In Ny vs. t at a relatively low T in (b) . The former display the deviation from 
the 'difussive' law. The latter confirms the expected logarithmic correction 
to a power law decay. The T-dependence and, especially, the saturation observed 
at high T are discussed in the text. 

A more precise determination of the growing length at temperatures close 
to Tkt should be given by the evolution of the density of free vortices only, 
pI. The identification of free vortices involves, however, some ambiguity. We 
used a simple-minded form that consists in counting as free all vortices and 
antivortices that do not have a neighbour of opposite vorticity at Euclidean 
distance r < Vc- By varying the parameter Vc from 1 to 2 we found that the 
second choice gives very good results, shown in Fig. din the form A^/lnA^/ 
against t for T = 0.4, 0.6, 0.8. (In [51 Vc is chosen to be equal to ^(t,T).) 
Consistently, ^ at sufficiently low temperatures (say, T < 0.4) but the 
two depart at higher temperatures, with the former yielding the correct out of 
equilibrium correlation length and being the one that ensures dynamic scaling. 



4.2.1 Temperature dependence 



The T-dependence in A and to in Eq. (12) has not been fully established yet. 
The assumption that at finite temperature the dynamics on scales that are 
shorter than ^(t,T) should be described, in the limit of large ^(t,T), by renor- 
malized spin- wave theory [33l |45] suggests A(T) = ps{T)V{T) and to{T) = 
a^/[ps(T)r(T)] at all T < Tkt- cl is the lattice spacing or microscopic cut-off, 
ps the renormalized spin-wave stiffness and F the renormalized kinetic coeffi- 
cient in the time-dependent Ginzburg-Landau model. The numerical data in [37] 
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Figure 3: (Color online.) (a) Time-dependence of the density of free vortices, 
expressed in the form A^/(t, T) In A^/(t, T) against after an infinitely rapid 
quench from Tq ^ oo to < T < Tkt given in the key. (b) Time-dependence 
of the growing correlation length, defined as = for values of T given 

in the key. Vortices and antivortices are considered to be free when the distance 
to their closest neighbour of the opposite vorticity is larger than a cut-off that 
in this figure has been chosen to be Tc = 2 lattice spacings. 



analyzed with an effective power- law C(t,T) A(T)t^/^(^) yield an effective T- 
dependent exponent that slowly increases with decreasing T. This is consistent 
with a microscopic time-scale to(T) in Eq. ([l2| that also weakly increases upon 
decreasing temperature. This in turn implies that A(T) is a weakly increasing 
function of T. We analyzed the T-dependence of the parameters A and by 
plotting 

i^{t,T)\iit against t (16) 

m Fig.||(b) and[4](a), using (^l and ^2 , respectively. Since the microscopic time 
to cannot be determined precisely due to the relatively long transient needed to 
reach the asymptotic law, and the fact that it yields a sub-dominant contribution 
anyhow, we fix it to one and determine the remaining free-parameter A. The 
data fall on straight lines with slope A(T). For all temperatures simulated, 
A(T) is an increasing function of temperature. As the temperature approaches 
Tkti ^{T) as obtained from the density of all vortices decreases, see Fig. [4j 
As already explained, this is an artifact of not having distinguished free from 
bound vortices. The agreement between the trend obtained from ^2 and is 
very satisfactory. The dashed line is a linear fit to the datapoints, A(T) Aq +cT 
with Ao — 1.2 and c = 6.8. 

The T-dependence of to and A are consistent with Berthier et aFs recognition 
that the Monte Carlo generated vortex configurations at a given waiting-time af- 
ter a quench to different temperatures are rather different. At high T (although 
below Tkt) the vortices tend to form pairs, as they should in equilibrium below 
TxT, while at lower T there are more free vortices. The suggestion is that the 
dynamics are faster at higher T (since vortex diffusion is faster) than at lower 
T 142J. 
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Figure 4: (Color online.) Temperature dependence of the parameter A(T) in the 
growth law (12) after an instantaneous quench from infinite temperature into 



the ordered phase, (a) The growing length determined as ^2{t^T). The slope in 
the long-time limit is A(T) and it increases with increasing T for all T's studied. 
The time-scale to{T) cannot be determined due to the long transient, (b) A(T) 
using four determinations: ^2, the total density of vortices and the density of 
free vortices with Tc = 1 and Tc = 2. The agreement between the first and the 
last one is very good. 



4.3 Spatial structure 

Further information on the vortex anti-vortex pair structure is obtained from 
the analysis of pair correlation functions that ignore or take into account the 
polarity of the vortices. In the former case we define a variable rii that equals 1 
if the site i is occupied by a defect irrespective of its polarity and otherwise. 
In the latter the variable Ui takes three values; = 1 if the site is occupied by 
a vortex, Ui = —1 if it is occupied by an anti- vortex, and = if there is no 
defect. The correlation function is then defined as G{r^t) = {ni{t)nj(t)) where 
r =\fi — fj\m. both cases. The first type of correlation does not yield additional 
information on the behaviour of the system with respect to C(r, t) studied above. 
Instead, the latter has a minimum at the preferred distance between vortices 
and anti- vortices |50]. In Fig. [s] we study G(r, t) at several instants after the 
quench, showing in (a) the raw data in the form rG(r, t) against r and in (b) 
the scaling plot (_l{t^T)rG{r^t) against r/(_y{t^T) with (^y{t^T) the correlation 
length obtained from the analysis of the density of free vortices (similar results 
are obtained using ^2)- The data confirm that this quantity is well described 
by dynamic scaling. The scaling function has a minimum at r 0.48 and 
it oscillates with very small amplitude around zero at long distances, a feature 
observed experimentally [49 but not captured by the approximate theory in [50 . 
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Figure 5: (Color online.) Polarity correlation function, (a) rG{r^t) against 
r at different times t after the quench into the low-temperature phase, (b) 
Scaling plot ^^{t)rG{r,t) against r/^y{t) with ^^(t) the length extracted from 
the analysis of the density of free vortices. The numerical data have been 
smoothened. 



4.4 Pair distribution function 

MC simulations give us access to the distances between vortices and anti- vortices 
and, in principle, to their pairing. The assignment of pairing is, however, a hard 
problem not free from ambiguity. We used the following simple algorithm. Given 
a configuration, we first computed and ordered all distances between vortices 
and anti- vortices. However, the lattice structure implies that, for sufficiently 
high density, some vortex (or anti- vortex) could be at equal distance from two 
(or more) anti- vortices (or vortices). In these cases we chose the pairing at 
random and we continued the procedure with the remaining defects. With 
this method we are not sure of finding the optimal pairing but, statistically, 
we expect all these pairings to be equivalent as the distribution of distances 
is concerned, a feature that we verified numerically. An example of the pair- 
distance construction is shown in Fig. [6j In panel (a) the configuration at 
t = 100 MCs after a quench from Tq ^ oo to T = 0.4 is represented as a 
Schlieren pattern in which a gray scale is proportional to sin^(20) with (p the 
angle formed by the local spin and the x axis. Each vortex emanates eight 
brushes of alternate black and white color. In panel (b) the vortices and anti- 
vortices are shown with (red) pluses and (blue) squares. A few pairings, as 
obtained with our simple algorithm, are highlighted in grey. 

Figure [7| displays the evolution of the distribution of pair distances after a 
quench from infinite temperature to T = 0.4. The log-log plot in panel (a) 
demonstrates that the initial probability distribution function (pdf) is a power 
law, with r ~ 3.5. Similarly to what was found in the study of geometric 

properties of coarsening in the Ising universality class [52 , small linear scales 
- as compared to - change significantly in time while larger scales do not. 
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Figure 6: (Color online.) Defect configuration at t = 100 MCs after quenching 
the sample to T = 0.4 in the ordered phase. The configurations are shown 
as Schlieren patterns over the full lattice {L = 100) on the left panel. The 
gray scale is proportional to sin^(2^) with (p the angle formed by the local spin 
and the x axis. On the right panel, the vortex configurations at time t = 100 
are shown. Vortices are represented by (red) pluses and antivortices by (blue) 
squares. The total number of vortices in the system is Ny{t = 100) = 80. We 
graphically show a few pairing of vortices and antivortices. 

This means that the tail of the distribution remains the initial power-law for 
r ^ (t) while the weight on smaller scales decreases significantly as a function 
of time. Panel (b) presents the scaling plot, (^^Np{r^t) against r/^^(t), in which 
one confirms that dynamic scaling is well satisfied. The crossover at r ^v(^) 
is clear in the figure. The kind of oscillating features at very short distances, 
r 3, might be due to the discrete lattice. 

The data in Fig. [t] is to be compared to results shown in Fig. 1 of Ref. [30 . 
The latter were obtained using the Fokker-Planck equation for the density of 
pairs of separation between r and r -\- dr proposed in [53 . This equation is 
tailored to hold at T < Tkt^ it does not take into account the effect of free vor- 
tices, and it does not capture the logarithmic correction to the diffusive growing 
length. The quench is done from equilibrium at T^t, where the distribution is 
the power law with r 4.70, to a very low T. The qualitative features of 
the two sets of data are the same. However, the quantitative properties are not. 
The initial conditions, and hence the initial distributions, are different. The 
time-dependent part has a peak in our case and it is totally fiat in [30 . Our 
pdf satisfies scaling with respect to the correct growing length. 
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Figure 7: (Color online.) Number of vort ex- ant i vortex pairs, A/p, as defined 
in the text, after the quench from T ^ oo to T = 0.4. (a) Np as a function 
of their distance r at different times given in the key. The dashed line is a 
guide-to-the-eye representing the power r~^-^. (b) Scaling plot. 

5 Dynamics after a slow cooling 

In this Section we discuss the evolution of the system after a slow cooling from 
high T to its low T phase. The cooling procedure used is graphically represented 
in Fig. [S] and it consists in the following. We suddenly quench the sample from 
T ^ oo (random initial condition) to Tq. In practice, we chose Tq = 2Tkt- 
We let the system evolve at this temperature a large number of MCs, typically 
trei = 2 X 10^ to 6 X 10^ depending on the system size. Equilibration occurs 
relatively rapidly far from Tkt so for all purposes we can assume the system 
to be in equilibrium at 2Tkt- We then cooled the sample by changing the 
temperature linearly in time according to 

T{t)=TKT{l-t/TQ) -TQ<t<TQ. (17) 

Within this time reference, annealing from 2Tkt starts at t = —tq. At negative 
times t the system is above the phase transition, at t = it reaches Tkt and 
at positive times it enters the ordered phase, tq is the inverse cooling rate. 
At t = Tq the environmental temperature vanishes and times are naturally 
bounded by tq. See Fig. [Sjfor a graphical representation of these cooling rate 
procedures. Although more complicated cooling procedures are also of interest, 
and have been analyzed in detail in the Id Glauber Ising chain [13^, we shall 
not discuss them here. 

5.1 Snapshots 

A first intuitive understanding of the defect density left over in the ordered phase 
after a slow annealing through the phase transition is obtained by looking at 
the spin configurations. 
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Figure 8: (Color online.) Sketch of the cooling procedures used. 



In Fig. [9] we display the system configuration on a section of linear size 
L = 100 of a system with size L = 256 at t = 200 MCs after having crossed the 
phase transition with two inverse cooling rates tq = 256, 2048 MCs. Note that 
the temperature at which the snapshots are taken is different, T{t = 200, tq = 
256) = 0.2 [(a) and (b)] and T{t = 200, tq = 2048) = 0.8 [(c) and (d)]. The 
figure demonstrates that the density of vortices increases with increasing tq, 
with Ny = 186 in the first case and Ny = 352 in the second. 



In Fig. 10 we show four configurations at t = rg, when the temperature has 
reached the value zero, with tq = 128, 512, 1024, 2048 MCs from (a) to (d), 
respectively. It is clear from the figure that the density of vortices decreases 
with increasing tq. We measured Ny = 180, 92, 68, 42 for increasing tq. In 
the following we shall make these statements quantitative. 

5.2 Equilibrium — out of equilibrium crossover 

The annealing procedure occurs by the following steps. In the high temperature 
phase the system is plagued with free vortices and ant i- vortices. The initial state 
thus contains a large number of free defects. Let us assume that the system is 
in equilibrium at Tq (= 2Tkt in our simulation). Under the cooling procedure, 
if this is slow enough, it evolves adiabatically, i.e., it remains instantaneously 
in equilibrium at T{t). The equilibrium correlation length is swept in time 
according to 

Ce,W^a«e''«l*/-«l"'' (18) 

and the density of vortices, estimated from py{t) £,eq{t)^ decreases exponen- 
tially: 

p^it) ^ af e-^'^Wr^r\ (19) 

However, this regime cannot last for ever, it crosses over to a different, much 
slower one, when the system falls out of equilibrium due to the fact that the 
relaxation time becomes much longer than the parameter variation imposed by 



the cooling rate procedure. Indeed, using Teq as derived in (15) and replacing 
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Figure 9: (Color online.) Comparison between the defect configuration at t = 
200MCs< min(rQ) after annealing the sample into the ordered phase with two 
different inverse cooling rates tq = 256 MCs and 2048 MCs. We show parts 
of linear size 100 of the full lattice L = 256. The configurations are shown as 
Schlieren patterns on the left [(a) and (c)]. The gray scale is proportional to 
sin^(20) with ^ the angle formed by the local spin and the x axis. On the right 
side, the vortex configurations at time t = 200 are shown [(b) and (d)]. Vortices 
are represented as (red) pluses and antivortices by (blue) squares. The total 
numbers of vortices in the system are Ny = 186 for the fast cooling in which 
T = 0.2 is reached and Ny = 352 for the slow one in which the configuration is 
measured at T = 0.8. 
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Figure 10: (Color online.) Comparison between the defect configuration at 
t = tq after annealing the sample to T = with four different inverse cooling 
rates tq = 256 (a), 512 (b), 1024 (c) and 2048 (d). We show a section of linear 
size 100 of the full lattice with L = 256. The total number of defects decreases 
with increasing tq and are Ny = 180, 92, 68, 42 from (a) to (d). 
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T{t) by its time-dependence, one finds 



zb^\t/TQ\- 



'\t/rQ\- 



(20) 



The system can no longer follow the pace of evolution set by the changing 
conditions. The crossover is estimated as the moment when the relaxation time 
reaches the characteristic time of variation of the temperatur^ AT/d^ AT ■ 



I eq 



—t ^1. The crossover occurs above the critical temperaure, AT = T — Tkt > 
and t > 0. The negative time —t is computed by assuming that the crossover 
occurs close to the transition, so that the critical scaling (l9| determines the 



equilibrium relaxation time Tgg. One then has | — t 
i :^ ar {i/rq) 



ji—t) and 



(21) 



The tq dependence of t can be obtained numerically through the Lambert W 
function, ie., the inverse of f{x) = xe^ . An analytic argument to get the leading 
Tq dependence relies on recasting (21) as 



(1 + i^) Inx + hi{TQ/ar) — b^zx 



with 



i/TQ . 



(22) 



Assuming that the solution is such that x ^ with rg/a^- ^ cx) we find, in first 
approximation. 



t 



Tkt \\n{rQ/ar) 



l/u 



(23) 



that is consistent with the assumption, t is smaller than tq - due to the log- 
arithmic correction - and this is in agreement with the fact that the system 
should fall out of equilibrium after the beginning of the annealing procedure at 
t = —Tq. Plugging now the 'first order' solution (23) in Eq. (22), we find a 
further logarithmic correction: 



t 



T-Tki 
Tkt 



b(:Z 



In 



TQ / h^Z \ 

ar \\n{TQ/ar) J 



(24) 



We compared the numerical solution to Eq. (21 ) to the estimate (24) for different 
Tq. For small and b^ \ the error is less than 10% for tq ^ 200. Taking 



into acount the full form (24) the correlation length at the crossover reads 





1 






_ln(rQ/a^)_ 





(25) 



-"^In the IdIM with Glauber dynamics, a model with an exponential divergence of the 
relaxation time close to the critical temperature Tc = 0, Krapivsky estimates —t for generic 
cooling rates by comparing the time needed to reach the critical point to the relaxation 
time [13) . For a linear cooling the two definitions are identical. 
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The length ^(t(T),T) coincides with the equihbrium correlation length for 
T(t) > f, when ^(£,f) reaches i, and ^(t(T),T) < | for T(t) < f. The 
crossover length, ^, increases and the crossover temperature, T, decreases with 
increasing tq. Consequently, one can get closer to the critical point evolving 
in equilibrium, and thus reach longer equilibrium correlation lengths, for longer 

For the sake of comparison, recall that in a conventional second-order phase 
transition with power-law divergencies, these quantities are both power laws [15 : 

i/TQ = (T - Tkt)/Tkt ^ (rQK)-'/^'+""-^ and i ^ a^(rgK)"/^'+""-^. i is 
also smaller than tq in this case. 



5.3 Numerical estimate of the crossover length 

In order to put the prediction for ^ to the test we compared the equilibrium 
and growing length extracted from the numerical data. We calculated the equi- 
librium values using systems with small linear sizes, L = 16, 20, 30, 50, that we 
let equilibrate. In Fig. [ll] (a) we display equilibrium data for ex pv with 
points. The data for T > 1.2 do not show finite size effects and are in agreement 
with the ones in [35] (obtained using much larger samples, L = 512) while at 
lower-T the spread of data is severe and the maximum values found for L = 50 
are way below the values estimated in [35 . The BKT exponential divergence is 
reported on the figure with solid line. The dynamic data for ^^[T(t)] obtained 
by using different cooling rates, tq = 128, 512, 2048 MCs, are shown with thin 
lines in both panels. Finally, the T- values at which the dynamic data should 
diverge from the static ones are signaled with dashed vertical lines. They are 
obtained from the numerical solution of eq. (21 ), by using the stemming from 



the fit of the BKT exponential divergence. The parameter a^- is, on the other 
hand, determined in such a way that the resulting T- values give a good estimate 
of the divergence of the dynamic data from the static ones. 



In Fig. 11 (b) the numerically determined t data are compared to the analytic 
estimate ( |24| ). The agreement between Monte Carlo data and analytic estimate 
is very good within our numerical accuracy. 



5.4 Out of equilibrium evolution 

At time —i close to criticality, the system falls out of equilibrium but it continues 
to evolve further, from T = T{-t) > Tkt to T(tq) = following the out of 
equilibrium critical dynamics rules. In the case of a finite rate quench, the total 
time spent in a critical region includes time t, which the system spent above 
Tkt after falling out of equilibrium. That is to say, the total time evolving with 
critical dynamics is At = t -\- 1. Then, we propose to match the equilibrium 
relaxation above T with the critical one below T with the asymptotic scaling of 
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Figure 11: (Color online.) (a) The equilibrium correlation length, at 
T > Tkt^ for different linear system sizes given in the keys are shown with 
points. The BKT exponential divergence, as obtained by fitting the data for 
temperatures higher than T = 1.2, where we do not have strong finite size ef- 
fects, are shown with think (red) lines. The time-dependent length ^^[T(t)] for 
different cooling rates indicated on the figure are shown with thin lines. The 
predicted Ts, for different tq, are signaled with dashed vertical lines. See the 
text for more details, (b) The crossover time t as a function of tq. The data 
points have been estimated from the numerical simulation - see panel (a) - and 
the line is the analytic prediction in eq. (24). 



the growing length 

_ r L,[T{t)] ^ t < -i{TQ) 

^^^^ {mt)] MAW)i) }'^^ ^ ^i-T(t) t > -i{TQ) 

Note that we assumed that the T-variation of the parameters A and to should 
be sufficiently slow so that we can simply replace their time-dependent values in 
the growth-law after a quench. We shall further fix to = 1 as in the quenching 
case. At At = 0: T{t) = T{—t) = T, £,eq{t) = C these formulae match. 

Let us confront the order of magnitude of the two terms in the second line of 



eq. (26) for times of the order of the inverse cooling rate, t ^ tq. Using Eq. (25), 
the first term is ^ :^ [tq/ (Irtq)^^^^^^^]^^^ . Since i is much smaller than tq, 
in the long tq limit At tq. The A prefactor is a function of T and hence of 
t but it should be bounded by the finite and non-vanishing limiting values in 
Fig. |4j Therefore, apart from the finite prefactor originating in A, the second 
term is of the order :^ [tq/ XivtqY^^ . Using v = 1/2, one concludes that the 
second term, describing the out of equilibrium dynamics below T, dominates in 
the large tq limit. (The same occurs in the conventional second order phase 
transitions although the mechanism is different [20j.) This fact can be observed 
in panel (a) of Fig. [12] where the growing under annealing with seven cooling 
rates given in the key is shown as a function of the time-varying temperature. 
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In Fig. 12 we present data for the dynamic growing length for different 
cooUng rates. Panel (a) displays the bare data for seven inverse cooling rates tq 
given in the key. Panel (b) is scaling plot that tests the hypothesis in Eq. (26) 
in the low-temperature region during the annealing process. We use the values 
of the parameter X[T{t)] obtained from the linear fit of data for the quench in 
Fig. |4j Before rescaling, the correlation length ^v{t) grows with time, i.e., with 
decrease of temperature, in such a way that at a given T, the value of increases 



with Tq (see panel (a)). After rescaling with the proposed asymptotic form (26) 
for t > —t{rQ), the values of ^ fall on top of each other. Perfect agreement with 
the analytic prediction would be the flat master curve ^v(^)/6owt(^) = 1- The 
master curve deviates from this prediction but the range of variation is relatively 
small, varying from 1 to 1.6, circa. We ascribe this mismatch to potential non- 
linear corrections to X[T{t)] and small variations in to[T{t)] that we have not 
taken into account. 
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Figure 12: (Color online.) Growing correlation length ^v{t) during annealing 
from 2Tkt to T = with different cooling rates tq. (a) Bare data for seven 
coolin g r ates given in the key. (b) Rescaling with the proposed asymptotic 
form (pel for t > -i{rQ). 



5.5 Numerical measurements of vortex density 



In Fig. [13] we study the density of all vortices during annealing. We present the 
data obtained by using five cooling rates in a log- log plot. The translation of 
the X axis by tq is done to give a common origin to all curves as well as to allow 
for a direct comparison with the results in [30 . 

Four regimes are evidenced in the figure. At very short times, as compared 
to a function of tq, the data remain close to the initial and rather large value 
py c:^ 0.2. This regime crosses over to one with a rapid relaxation in which 
the density of vortices decreases significantly (the shoulder in the curves). The 
crossover occurs dit t ^ —i. The relaxation continues in a smoother way and 
it is close to a power-law decay with logarithmic corrections as expected. This 
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regime ends when the curves reach the hmit of the instantaneously quenched 
data to T = 0.4 that they subsequently follow until t = tq and T = 0. In our 
setting the density of defects increases with increasing tq at fixed t-\-rQ. 

The qualitative features of our data are very close to the ones obtained by 
Chu and Williams [30 for a quench or annealing from equilibrium at Tkt to 
below the critical point. As already said, their method is to solve numerically 
the Fokker-Planck equation for vortex-pair dynamics in conjunction with the 
Kosterlitz-Thouless recursion relations derived in [53 . This approach can only 
be used for initial conditions at and below Tkt and does not include the effect of 
free- vortices. In consequence, the scalings found are consistent with power-laws 
without logarithmic corrections. 
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Figure 13: (Color online.) Time-dependence of the total number of vortices 
after annealing to T = with different inverse cooling rates tq given in the key. 
The data are shown with points in the form Ny \nNy as a function of t -\- tq. 
With a (red) line we plot the number of vortices after an infinitely rapid quench 
to T = 0.4 < Tkt- 



Finally, we extracted the remanent number of vortices at the end of the 
cooling procedure, t = tq and T = 0, and we plot it as a function of tq in 
Fig. 14 The data correspond to the direct counting of (free) defects. The total 
number of defects varies from 10^ to 10 circa. On the analytic side, the scahng 



assumption in Eq. (26) suggest 



Pv{tq) 



'Q 



'Q 



V ^ in(At(rQ)/to[T(rQ)]) 



r 



+ A 



A 



ln[(TQ 
{rq + 1) 



1 



(27) 



ln[(TQ+t)Ao] 

where A and to in the last expression are the zero-temperature values. We 
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M] by using 
and we found good agreement 



compared the above assumption with the n ume rical data in Fig 
the values of i obtained numerically in Sect, 
within our numerical accuracy. The fact that the data points bend at short 
tq is well captured by the addition of t in the numerator. This term becomes 
negligible for longer tq. At intermediate values of tq the decay can be easily 
confused with an effective power-law, t'^'^^ ^ shown in the figure with a dashed 
(blue) line. For longer values of tq the data confirm the logarithmic correction. 
Note that the decay is slower than in the discrete symmetry breaking case with a 
low-T phase with trully long-range order, where one finds py{t = tq) c:^ Tq^ [20^ . 
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Figure 14: (Color online.) rg-dependence of the density of vortices after anneal- 
ing the system to T = 0. Points represent the numerical data obtained from the 
direct counting of topological defects, the dashed (blue) line shows an effective 
power law Tq^'^^, while the full (black) line corresponds to the prediction in 
Eq. (27). Similar results are obtained from py (not shown). 



6 Conclusions 

This paper should contribute to the understanding of the dynamics of topolog- 
ical defects across classical thermal phase transitions. 

The framework we focused on is one in which the system of interest is coupled 
to an environment in equilibrium at a given temperature. The contact to the 
bath induces dissipative stochastic dynamics that, for concreteness, we chose not 
to conserve the order parameter. The system is taken across a phase transition 
by either varying the temperature of the external bath or by tunning a parameter 
in its Hamiltonian. The disordered (high temperature) phase is plagued with 
topological defects. After entering the ordered (low temperature) phase the 
system's configuration is out of equilibrium and the system relaxes by changing 
their organization and reducing their density. The nature and structure of 
topological defects depends strongly on the dimension of the order parameter, n. 



24 



and the dimension of space, d. For the scalar case, n = 1, the defects are domain 
wahs and these are points in d = 1, hnes in d = 2 and surfaces in d = 3. For 
n = d the defects are point-hke and the example at hand are vortices in d = 2. 
While the ordering kinetics of systems subjected to instantaneous quenches has 
been studied at length with field theoretical and numerical techniques [9" , the 
analysis of the effect of infinitely slow cooling rates has not been developed to 
the same extent (see [lOl [11] [121 ^M)- The Kibble-Zurek mechanism - based 
on critical scaling above the critical point - provides a concrete prediction for 
the density of defects left over across a phase transition when this is crossed 
at very low rate. This proposal captures correctly the moment and parameter 
value at which the system falls out of equilibrium (see p.3 for a very detailed 
calculation in the Id Glauber Ising chain) but neglects, incorrectly in dissipative 
thermal phase transition, the relaxation dynamics in the ordered phase. In [20J 
the case of a second-order phase transition with discrete spontaneously broken 
symmetry was discussed. In this paper we studied a different type of transition, 
the BKT one in the 2d XY model. 

To summarize our results, we first analyzed in more detail than previously 
done the out of equilibrium relaxation of a 2d XY model infinitely rapidly 
quenched from Tq ^ oo to T < Tkt- Using MC simulations we confirmed the 
functional form of the growth law ^(T,t) [A(T)t/ Int/to]^/^ derived in [32, 33] 
and we determined its temperature dependence numerically. Although the log- 
correction has a theoretic foundation, this law has been frequently confused with 
an effective power law, especially in the field of freely decaying 2d turbulence as 
modelled with a Ginzburg-Landau approach [4j. We also analyzed the evolu- 
tion of the structural properties such as the distribution of pair distances or the 
polarity correlation function in the course of time and we found that they all 
satisfy dynamic scaling with respect to the same growing length. Importantly 
enough, we demonstrated that the out of equilibrium dynamic properties are 
determined by the density of free vortices in the low-T phase as opposed to 
the density of all vortices that includes the equilibrium contribution of bound 
pairs. This fact is the reason for the failure of dynamic scaling when the length 
determined from the decay of the total density of vortices was used [46* i47j • 
Next, we studied cooling rate dependencies, we proposed a scaling form for the 
growing length under these circumstances, and we checked it numerically. We 
found that the density of free topological defects not only depends on the cool- 
ing rate used but does also on the time spent in the low temperature phase, 
where vortices and anti-vortices tend to bind or annihilate. In particular, the 
density of vortices at t :^ rg is py [rg/lnrg]"^ for very long tq. In short, 
the dynamics in the critical low-T phase of the 2d XY is crucial to determine 
the density of topological defects. 
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